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1. Introduction 


There are a wide variety of computational tools available for analysis and design of munitions 
and armor. While simulations of high-rate dynamic loading have predominantly used finite 
element and hydro codes, emerging methods also receive attention. One relatively new approach 
is peridynamics (7-3). Initially proposed at Sandia National Laboratories (7), it was developed 
with fracture simulations in mind. The concept is that integrating over pair-wise interactions 
among material points ameliorates some of the difficulties encountered when fonnulating 
fracture problems with continuum-based methods that rely on continuous fields to evaluate 
gradients. Most of the use of peridynamics to date has been for studies of fracture and 
fragmentation ( 4 - 6 ). As the peridynamics codes now include plasticity ( 2 , 7) and are becoming 
more widely available ( 8 , 9 ), there is a desire to apply them quantitatively and to a broader class 
of problems. 

As with many tools, there will be some areas where peridynamics will be advantageous, and 
there may be other areas where it is not the most appropriate method. The intent of this report is 
to begin defining the role for peridynamics by providing an assessment for high-rate defonnation 
simulations, including large strain deformations that often precede facture and material flow 
accompanying ballistic impact. 

The report will begin with an overview of the peridynamics equations; first the micro-elastic and 
micro-plastic models will be outlined, and then the newer state-based formulation will be 
presented along with a comparison to traditional finite elements. Numerical examples will be 
presented to assess peridynamics results in deformation scenarios that typically accompany 
fracture or ballistic impact. The origins of any unexpected behaviors are explored through 
parameter studies and analysis. Recommendations will be made regarding the range of 
applicability of peridynamics for high-rate deformation problems. 


2. Peridynamics Formulation 


The initial fonnulation of peridynamics (7) was focused on integration of pair-wise particle 
interactions over some domain surrounding a material point. This provides a convenient 
framework for elastic models. However, the pair-wise model construct is not sufficiently flexible 
to incorporate traditional nonlinear constitutive relations, nor is it amenable to representing field 
relations that are typically expressed as differential equations. The peridynamics framework was 
subsequently extended to a state-based approach ( 2 , 7 ) to facilitate use of common nonlinear 
constitutive models and more general field equations. Both the micro-elastic and the state-based 
formulations are discussed. 
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2.1 Micro-elastic and Micro-plastic Models 

Peridynamics solves the traditional momentum equation at material points within a body. 

f(x, t) + b(x, t) = pii (1) 

Here t is time, x are the coordinates of the point, u is the displacement of the point, / is the 
force acting on the point by the surroundings, b is the body force, and p is the density. The 
double dot denotes the second time derivative. The force is obtained from a function of pair-wise 
differences in displacements and coordinates integrated over the region 91 surrounding the point. 

f(x,t)= J 3[u(x',t) — u(x,t),x'— x]dV x , (2) 

'tv 

The prime denotes the other points in the region connected to the point of interest. To make the 
solution tractable, the integration over material points within region 91 is replaced by a 
summation over more widely space nodes within a radius called the horizon. The horizon is set 
to 3.015 times the grid spacing by default in the peridynamics codes EMU (5) and KRAKEN (9). 

Failure occurs by progressive fracturing of individual li nk s as the local strain limit is exceeded. 
The failure criterion is similar to that available in many finite element codes, but the progressive 
failure aspect is less common. Once a link has failed, the connection between the nodes is lost. 

The micro-elastic and micro-plastic models both provide a linear force-displacement relationship 
between nodes, with a cap on the peak force (yield) and a critical force to failure. The response 
in unloading is different for the two treatments in that the elastic model unloads along the 
original path and the plastic model unloads with the elastic slope. 

The micro-elastic/micro-plastic model is limiting in two respects. There is a single elastic 
constant, whereas an isotropic elastic continuum material has two. For the micro-elastic and 
micro-plastic models, the modulus is related to the speed of a longitudinal stress wave and the 
Poisson’s ratio is fixed at 1/4 (7). The Poisson’s ratio cannot be varied. The second limitation is 
that yield applies to the total model behavior. In a traditional elastic-plastic continuum model, the 
yield affects the deviatoric part of the behavior while the volumetric part remains elastic. 

2.2 State-based Model 

The state-based peridynamics fonnulation addresses some of the limitations of the micro-elastic 
and micro-plastic models. It is useful to begin the state-based peridynamics presentation by 
giving some of the basic equations in comparison to Lagrangian finite element methods 
(FE/FEM). Operationally, both peridynamics and Lagrangian FEM have nodes attached to 
material points that track displacements and velocities. In peridynamics, the strain and stress are 
also resident on the nodes, and quantities are obtained by integration over other nodes within a 
given radius called the horizon. There is no explicit mesh. In FEM, elements are constructed with 


2 



the nodes as the vertices, and the strain and stress are obtained by integration over the element 
volumes. 

2.2.1 Kinematics 

The expressions for the continuum representation of the average deformation gradient in 
peri dynamics over its horizon (7) and the average deformation gradient within a finite element 
(10) are, respectively 

¥ Peri = J <u(|f|)Cr 0 0 dV x < K- 1 and F fe = ±J (^) dV . (3) 

Both methods involve volume integrals, x are the coordinates of a material point initially located 
at X. The relative position between two material points in the initial configuration is define as 
^ = X' — X, and the relative position in the current configuration is Y = x' — x. The volume 
integration in peridynamics is over all points within the horizon except the focus point, x. The 
weight, <y(|f|), is generally equal to one except when bonds are fractured, and the shape tensor 
is defined as 


K= [ <o(\mf®OdV x ' . (4) 

Heuristically, for a unifonn point spacing in peridynamics, the integrated volumes from 
equations 3 and 4 cancel, as does a ^ from each. This leaves the sum of the current distances 
between points divided by the sum of the original distances. This result has the appearance of the 
definition of a derivative, as it must to produce the defonnation gradient. Hence, in a loose sense, 
the defonnation gradient in peridynamics is produced by the definition of a derivative while it is 
a fonnal mathematical derivative in the finite element formulation. A crack in the continuum FE 
formalization is a singularity, while in peridynamics it is a spike that can be made arbitrarily 
large or small depending on the chosen material point spacing. 

This distinction in the integrand singularity becomes academic when considering the numerical 
implementation of the techniques on a discrete grid of nodes. In discrete form, the deformation 
gradient for the two methods from equation 3 becomes 

K N 

F Peri = 2[<y(|f„|)(K„ 0 f n )V n ] K 1 F fe = ^ n ®B n . ( 5 ) 

n=1 n=1 

The B n in the FE representation, equation 5b, is the standard difference operator. The summation 
for peridynamics is over the nodes in the horizon, and for finite elements is it over the nodes 
defining the element. The discrete relative position vector is = X n — X m , where X m is the 
initial location of the focus node, and the definitions of Y and K follow similarly. Recognizing 
that K -1 is a constant in equation 5a, it can be taken inside the summation and rewritten as 
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(6) 


K K 

F p eri = Y j [Y n ® (.fn-K-'aQfnWn)] = £ [(*„ ~ X m ) ® B n ] , 

n=1 n=1 

where the B n is a vector of dimension inverse length operating on differences between 
coordinates. 

As an illustration of what is calculated numerically, consider the 2D deformation gradient 
computed for the center point of a set of nine regularly spaced points, figure 1, with both x- 
spacing and v-spacing equal to d. For peridynamics the horizon is set at 1.5, which is smaller 
than typically used, but it facilitates comparison with the FE. The defonnation gradient from the 
FE solution is averaged for the four elements surrounding the center node, which is also not 
normally done, but it likewise facilitates the comparison. 



Figure 1. Node numbering for example in 
equations 7 and 8. 


The 2D deformation gradient components computed for the two methods are 

Xi + x 3 — x 4 + x 5 — x 6 + x 7 ) (— x 1 — x 2 — x 3 + x 6 + x 7 + x 8 ) 

' 6 d L(-yi +y 3 -y4 + ys-y6 + y 7 ) (-yi - y 2 - ys + y 6 + y 7 + y 8 ) 

1 [(—X! + x 3 — 2x 4 + 2x 5 — x 6 + x 7 ) (—x 4 — 2x 2 — x 3 + x 6 + 2x 7 + x 8 ) 

8 d L(— y 4 + y 3 - 2y 4 + 2y 5 - y 6 + y 7 ) (- 


c ’Peri _ 1 [( 

r ij 


nFE 


'Xi 

■yi 


y 3 + y 6 + 2y 7 + y 8 ) 


( 7 ) 

( 8 ) 


In this example, the distinction between the two methods is in the weighting of points at the 
center of the edges of the grid. The operations on the coordinates are otherwise similar in this 
example. 


The number or nodes involved in the summations for peridynamics and FE simulations will be 
quite different in typical applications, potentially leading to significantly different computational 
results. The horizon for peridynamics simulations will be larger (the default is 3.015 times the 
average node spacing), and the weighted sums and differences in 3D could span more than 100 
nodes. This makes peridynamics a nonlocal method which has additional consequences for 
solving the fundamental equations. On the other hand, deformation gradients for unifonn strain, 
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hexahedral finite elements will be computed by weighted sums and differences of eight nodes. It 
is a local method. The similarity is that both methods use weighted sums of nodal values; they 
will not necessarily give similar results. 

2.2.2 Weighting of Distant Nodes 

Peridynamics is a nonlocal method in that the behavior at a material point can be affected by the 
deformation and stress fields beyond its local vicinity. The weighting function, u>, can be used to 
control the influence of distance points relative to points in the immediate neighborhood. 
However, the current defaults in EMU are to have the weighting factor equal to one for all points 
within the horizon prior to failure. The user’s manual (5) does not indicate a way to change the 
weighting at this time. 

With a weighting factor of one, the structure of peridynamics places a greater weight on nodes 
near the horizon than on nodes near the focus point. As an example, consider a unifonn 1-D grid 
with spacing, d, subjected to a uniform stretch, A. The horizon is 3.015 times d. This has initial 
relative position vector magnitudes, equal to d for the nodes adjacent to the focus node and 
equal to 3d for nodes near the horizon. The current relative position magnitudes are Ad for the 
adjacent nodes and 3 Ad for the nodes near the horizon. Multiplying these as indicated in 
equation 3 a results in the nodes near the horizon having nine times the weight of the nodes 
adjacent to the focus point when computing the sums. 

2.2.3 Stress and Nodal Forces 

State-based peridynamics and FE both use the same momentum equation, equation 1, and similar 
constitutive models to integrate Cauchy stress given the applied defonnation history. Formulas 
translating the stress to a peridynamics nodal force vector state (7)* and to nodal forces in FE 
(10) are, respectively 

KU Peri = <u(|f n | >m ■ K- 1 • SnVm and f F n E = B n -aV . (9) 

Similar to equation 6, and using the symmetry of K, equation 9a can be written in tenns of B n as 

ftfn) Peri = fn- K-'-amOiMnDVm = B n CT m . (10) 

As with the expressions for the deformation gradients, the relations for the forces have 
similarities. However, there are notable differences in how the individual forces are combined to 
give the net force on a node. Peridynamics uses a sum of differences between force vector states 
integrated over the volume. 

K 

f(Xm ) = ^[/On - *m) - f(x m - X n )]V n (11) 

n =1 


The force vector state defined in equation 25 of reference 7 appears to be dimensionally incorrect for a force since K has 
dimension of L~ 5 . The presumed missing volume has been included in equation 9 above. 
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The resulting nodal force involves a difference in stresses among neighboring nodes. These 
stress differences can be viewed as comprising a gradient in the stress field, which resembles the 
strong form of the momentum equations. Finite element formulations use a weak form of the 
momentum equations, and the forces from equation 9b are simply summed at the nodes. 

2.2.4 Failure Modeling 

Failure for the state-based materials in the peridynamics code EMU is triggered when a critical 
bond strain is attained. The bonds between nodes are broken individually as they reach the 
failure criterion, resulting in a progressive material failure. FE methods are similar in that 
connections among the nodes of an element are modified when the element reaches a failure 
criterion. The connections within an element can be weakened gradually or removed abruptly 
depending on the failure model used. The gradual failure progression in peridynamics would 
make fracture propagation less sensitive to numerical discretization than abrupt element failure 
in FE methods. 

After failure, the peridynamic bonds are broken. In tension, the nodes are free to separate without 
constraint. To prevent unrealistic densification of failed materials in compression, short-range 
interaction forces are imposed when nodes are closer than some critical distance. 


3. Numerical Simulation Results 


Peridynamics codes EMU (5) and KRAKEN (9) were obtained from Sandia National 
Laboratories and installed on the Department of Defense (DoD) high-performance computing 
platform (HPC), harold. Although the parallel executables were built, the simulations were 
sufficiently small enough that all runs were on a single processor. 

The simulations were designed to assess basic dynamic material response that would occur prior 
to or in combination with fracture. The goal was to examine quantitative aspects of physics that 
represent basic behavior, rather than features of the solution that would be highly variable and 
model dependent. Hence, fracture behavior itself is not considered. The simulations were as 
follows: low velocity ID wave propagation, inhomogeneous simple shear, and low speed impact 
with a spherical projectile. Since neither code incorporated artificial viscosities necessary for 
shock simulations, the loading rates were relatively slow, and oscillatory results are expected. 

The use of wave propagation examples is not to focus on wave mechanics, but rather to impose 
sharp gradients in the defonnation field in a configuration where the response is well-understood. 
There are many more general configurations that would result in steep gradients, such as 
plugging associated with penetration; but the response in these situations is not understood as 
well, and errors would be more difficult to quantify. 
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3.1 Micro-elastic and Micro-plastic Models 

The micro-plastic model is the same as the micro-elastic model except for the unloading 
response, so only the micro-elastic model was evaluated. The material used for the ID micro- 

■j 

elastic simulations resembles aluminum. The density was 0.0027 g/mm and the longitudinal 
sound speed was 5750 mm/ms. The yield strength and failure strain were varied as part of the 
assessment and will be specified for each simulation. 

3.1.1 ID Elastic Wave Propagation 

Elastic wave propagation was investigated using a body 4 mm x 4 mm x 1 mm thick and with a 
regular node spacing 0.05 mm. The horizon was 0.1575 mm, 3.015 times the grid spacing. The 
yield stress was 1 GPa to inhibit yielding, and the fracture strain was set to 1000 to preclude 
fracture. 

Both longitudinal and shear waves were simulated using KRAKEN (9). A 10 mm/ms 
longitudinal or transverse velocity was imposed on one 4 mm x 4 mm face at t=0, and the 
opposite face was stress free. Velocities were recorded at nodes at the surface where the velocity 
boundary condition was imposed; at depths of 0.30 mm, 0.35 mm, 0.40 mm, 0.45 mm, and 
0.50 mm; and at the free surface. The velocity histories are shown in figure 2a for the 
longitudinal wave and in figure 2b for the transverse wave. The waves are smooth and well- 
behaved. The oscillations are expected given the spatial discretization and no artificial viscosity 
to treat the shock dissipation. 

The deformation and forces are smoothed over the nodes within the horizon resulting in a 
spreading of the wave front. This smoothing was examined by rerunning the longitudinal wave 
simulation with a horizon of 0.071 mm, 1.42 times the grid spacing. The resulting velocity 
profiles are shown in figure 3a. Comparing this with figure 2a, it is noted that narrower horizon 
results in a sharper rise time for the wave front. The nonlocal nature of the peridynamics 
formulation is responsible for spreading the gradient. 

The wave propagation was analyzed quantitatively to obtain measures of the wave speed and the 
front width. Given the broad wave front width, these quantities will vary based on how they are 
measured. Using the data in figures 2a and 3a, times are determined when the nodes reach 10% 
and 90% of the imposed velocity. These data from the five internal nodes are plotted in figure 
3b. The free surface nodes were not included because the reflection from the free surface might 
modify the timings. The slope of the line connecting the two end points of each line is given; this 
is the wave speed at the particular location. 
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Figure 2. Longitudinal (a) and transverse (b) particle velocities at the indicated positions for the micro-elastic 
material with a horizon of 0.1574 mm, 3.015 times the grid spacing. Results for (a) and (b) are for 
distinct calculations where a longitudinal or transverse velocity, respectively, was imposed on the 
surface. 
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Figure 3. Longitudinal particle velocities (a) at the indicated positions for the micro-elastic material with a 
horizon of 0.071 mm, 1.42 times the grid spacing. Times (b) at which nodes achieved velocities of 
1 mm/ms and 9 mm/ms for the horizon values indicated in the legend. The calculated wave speeds are 
given in the legend. 

Several observations can be made from figure 3b. First, the wave speed is significantly greater 
than the value placed in the input deck. Second, for both horizon values, the wave speed at the 
90% velocity is less than at the 10% velocity. This implies that the wave is broadening as it 
traverses the sample. Third, using the measured wave speed, the average distance the wave 
traveled between 10% and 90% velocities is approximately 0.19 mm for the 0.1575 mm horizon 
and approximately 0.11 mm for the 0.071 mm horizon. These distances indicate the wave front 
width, which corresponds to roughly four and two node separations, respectively. 
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3.1.2 ID Elastic Wave Propagation on a Sheared Grid 

The elastic wave propagation simulations with the 0.1575 mm horizon were repeated using a 
grid where each successive plane of nodes in the thickness direction (1 mm direction) was 
shifted parallel to this plane by 20% of the node spacing. When sectioned, the nodes are the 
vertices of parallelograms, as shown in figure 4. 



Figure 4. Sheared grid used to assess the volume-shear 
coupling response indicated in figure 5b. 

The results of these sheared grid simulations are shown in figure 5. The longitudinal particle 
velocities, figure 5a, are similar to that of figure 2a, but the wave arrived at the free surface 
approximately 0.008 ps sooner. Figure 5b shows the node velocities perpendicular to the 
direction of wave propagation. The peak transverse node velocity is about 2% of the applied 
velocity. This volume-shear coupling is an anomalous result related to the node distribution 
within the horizon. If the grid shift is 0 or 0.5 of the node spacing, the configuration is 
symmetric, and the transverse velocity is zero. This shear will be explored further in section 
4.1.1. 
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Figure 5. Longitudinal (a) and transverse (b) particle velocities from a single ID longitudinal wave calculation 
with a sheared grid using micro-elastic material. The transverse velocity should be zero. 
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3.1.3 ID Plastic Wave Propagation 

The wave propagation simulations of section 3.1.1 were repeated using KRAKEN, with the yield 
strength lowered to 80 MPa to exercise the cap on the micro-elastic model. All other aspects of 
the input were identical to those in section 3.1.1. 

The resulting node velocity history plots for the applied longitudinal and transverse velocity 
boundary conditions are shown in figures 6a and 6b, respectively. The nodes interior to the 
specimen reach a velocity of approximately 4 mm/ms, with little oscillation due to the dissipative 
nature of the plastic deformation. The velocity of the interior nodes is governed by the stress cap, 
which limits the driving force for nodal accelerations. 

The critical point to note is that the interior velocities only reach 4 mm/ms velocity, while the 
surface is still moving at 10 mm/ms. In the longitudinal wave simulation, this indicates a 
continued densification over time. This result is not physically realistic, as the compressive 
response after yield should be governed by the equation of state. For the shear simulation, the 
difference between the surface and internal velocities corresponds to strain localization at the 
surface. 
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Figure 6. Longitudinal (a) and transverse (b) particle velocities at the indicated positions for the micro-elastic 
material with a stress cap of 80 MPa. Results for (a) and (b) are for distinct calculations where a 
longitudinal or transverse velocity, respectively, was imposed on the surface. The persistently lower 
velocity of the interior nodes in a) indicates a continuous volume decrease. 

3.1.4 Post impact Response 

The behavior of the micro-elastic model following compressive failure was examined by 
simulating a sphere impacting a brittle block. A 30 mm square target with a 0.5-mm grid-spacing 
was impacted by a dense, rigid, 10 mm diameter sphere traveling at 200 mm/ms. The target had a 
density of 0.0022 g/mm , a sounds speed of 5750 mm/ms, a yield strength of 1300 MPa, and a 
failure strain of 0.001. The initial configuration is shown in figure 7a. 
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The configuration 0.012 ms after impact is shown in figure 7b, where the false-color image is of 
material dilatation. The densification under the projectile at this time is approximately 40%, 
much greater than would be expected under these impact conditions. The links connecting the 
nodes are removed after failure, and the material is compressed. It is not possible to enforce 
volume constraints from an equation of state separately from material failure using the simple 
micro-elastic model, and the post-impact densification is regulated only by the short-range 
interaction forces. To produce an acceptable response in compression, the short-range 
interactions would have to be modified to correspond to an equation of state. This type of post¬ 
fracture response is not an issue in most tensile failures where the fractured interfaces do not 
interact. 



Figure 7. Dense, rigid sphere impacting a brittle material: (a) initial configuration and (b) dilatation after 
0.012 ms. 


3.2 State-based Constitutive Model 

The state-based models remove many of the limitations of the micro-elastic and micro-plastic 
models by providing a connection to traditional continuum material constitutive models. Several 
of the assessment simulations presented in the previous section are repeated with the state-based 
models in EMU. The properties for the wave problems are again chosen to resemble aluminum 

•5 

with a density of 0.0027g/mm , Young’s modulus of 70 GPa, and Poisson’s ratio of 0.3. The 
yield strength was set arbitrarily to provide either elastic or plastic defonnation. 

3.2.1 ID Elastic Wave Propagation 

Elastic ID wave propagation analyses on 1 mm thick specimens were run in EMU using a 1 GPa 
yield strength and imposing a 10 mm/ms velocity. The response with the state-based models is 
sensitive to the averaging within the horizon near the boundary (77), so the velocity boundary 
condition was applied to all nodes within the horizon of the surface. 
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Plots of longitudinal particle velocity at several locations through the sample thickness are 
shown in figure 8a. It can be seen that the wave propagation is not steady through the thic kn ess. 

The developers recommended (11) adding drag terms to control zero-energy modes and reduce 
the noise. The drag algorithm replaces the velocity of a node with a weighted sum of the original 
velocity and the average velocity of all the nodes in the family. Since there is no damage in these 
problems, the adjusted nodal velocity is given by 

K 

K = a-D)V n + D-Y j V m> ( 12 ) 

m =1 

where the drag coefficient, D, has been set to 0.05. Application of the drag algorithm reduced the 
variability in nodal velocity and produced a reasonably steady velocity profile for the nodes 
beyond the center (figure 8b). However, the drag algorithm is a smoothing technique and does 
not specifically target hourglass models. It also alters aspects of the solution not related to 
hourglassing. 



Figure 8. Longitudinal particle velocities obtained from an elastic state-based model (a) without and (b) with 
viscous drag and hourglass control. 

It is notable that the nodes at 0.9 mm and the free surface have nearly the same velocity-time 
history. Since this is the same with and without the drag algorithm, it is likely due to a 
combination of the averaging within the horizon, a nonlocal effect, and the long time period for 
diffuse wave reflections from the free surface. 

Wave velocity is more difficult to discern than in section 3.1.1, but it is approximately 
8000 mm/ms when the particle velocity is 0.1 mm/ms, and approximately 5400 mm/ms when the 
particle velocity is 0.9 mm/ms. The latter roughly corresponds to the expected longitudinal wave 
velocity for the material parameters provided. As with the micro-elastic model, the higher wave 
speed at the lower particle velocity indicates that the wave is spreading over time. 
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The drag algorithm in equation 12 averages velocities on every time step, resulting in significant 
smoothing of the solution. Since it is applied once per time step, the same calculation run with a 
smaller time step would have many more smoothing operations applied. This time step effect is 
illustrated by comparing figure 8b, which was run in 23 time steps, to the same input run with 
179 time steps, shown in figure 9. The results from figure 8b are included as thin dashed lines in 
figure 9. It is evident that the oscillations are damped further, and that the rise time for the waves 
has been lengthened by increasing the number of time steps. 



Figure 9. Calculation shown in figure 8b repeated with 
a smaller time step to illustrate the dissipative 
influence of the drag algorithm. The thin dashed 
lines are the curves from figure 8b. 

3.2.2 ID Plastic Wave Propagation 

Plastic ID wave propagation analyses on 1-mm thick specimens were run using the state-based 
model in EMU with a 50-MPa yield strength and imposing a 10-mm/ms velocity. As with the 
elastic solution, the velocity boundary conditions were imposed on several layers of nodes near 
the surface, and drag and hourglass algorithms were used to smooth the velocity fields. 

Particle velocities at several locations through the sample thic kn ess are shown in figure 10. 
Compared to the elastic solution, the rise time for the wave is considerably longer. The ki nk 
evident in the curves 0.5 mm and deeper is associated with the elastic-plastic transition. 
However, the response is smoothed considerably since the velocities are averaged over the 
horizon, and the stresses are effectively averaged over two horizon distances surrounding a node. 
This is in addition to the smoothing caused by the drag algorithm. The node spacing would have 
to be significantly finer to capture the nearly vertical rise and plateau normally associated with 
the Hugoniot elastic limit. 
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Figure 10. Longitudinal particle velocities for the 

state-based model where the yield strength 
is set low to allow plastic deformation. 


3.2.3 Sheared Grid 

The state-based model with a 50-MPa yield strength was run with a 10-mm/ms impact velocity 
using the sheared grid in figure 4 to determine if the volume-shear coupling was also present 
with the state-based models. Results are shown in figure 11 for simulations without and with the 
drag algorithm. It is evident that the incongruous volume-shear coupling is still present, but the 
drag algorithm imposes sufficient constraint to reduce the shear error by an order of magnitude. 



Figure 11. Transverse particle velocity for simulations using a sheared grid. Results are shown for calculation 
(a) without and (b) with the drag algorithm. 


3.2.4 Further Coupling Effects with Nonuniform Grid 

The effects of grid nonuniformity were further examined by propagating a wave through a 
spatially graded grid. A section through the center plane of the grid is depicted in figure 12. The 
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grid spacing is a uniform 0.05 mm in the vertical direction and in the depth direction. In the 
transverse direction, the mesh is unifonnly spaced at 0.0294 mm for 1 mm on the left, followed 
by 2 mm of graded mesh where the node spacing increases by 2.5% for successive columns of 
nodes. The last 1 mm is also a uniform grid with a spacing of 0.0833 mm. Wave propagation and 
shear calculations were run with this grid. 
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Figure 12. Simulated wave propagation through a spatially graded mesh. The plots show (a) longitudinal 
velocity (m/s), (b) transverse velocity (m/s), (c) longitudinal-transverse shear stress (MPa), and 
(d) longitudinal stress (MPa). 


The first simulation is a 50-mm/ms velocity imposed on the upper surface, resulting in a wave 
propagating vertically downward through the mesh. Symmetry boundary conditions are applied 
over 0.15 mm strips of material in the out of plane and in the lateral directions to suppress edge 
effects. Results are shown in figure 12. The longitudinal velocity and stress—figures 12a and 
12d, respectively—show the location of the wave front. They also show that the solution is 
significantly different in the center region of the specimen, where the grid spacing is graded, than 
on either end where the grid spacing is uniform. 
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A transverse velocity is induced by the mesh gradation, shown in figure 12b. The magnitude of 
the transverse velocity is roughly 15% of the applied velocity, which is a very large volume- 
shear coupling effect. This introduces a correspondingly large shear stress, which is captured in 
figure 12c. The transverse velocity and shear stress are spurious results related to the node 
distribution within the horizon. 

An additional grid effect evident in these simulations is the influence of mesh spacing orthogonal 
to the wave propagation. The velocity and stress are approximately 5-10% higher in the fine grid 
region on the left than in the coarse mesh region on the right. The grid is uniform on the left and 
right of the body, so this lateral node spacing effect is independent of the grid gradient effects. 

The second calculation with the same graded mesh examines the potential complementary shear- 
volume coupling. In this simulation, all of the nodes are constrained in the vertical and out-of¬ 
plane directions. A 50-mm/ms transverse velocity boundary condition is imposed on the bottom 
and top 1 mm of material. The center 1 mm of material is given a -50-mm/ms transverse 
velocity. This creates two 1-inm thick regions subjected to simple shear boundary conditions 
similar to plugging. The imposed x-direction velocity field is shown in figure 13. 

The resulting out-of-plane stress and dilatation are shown in figures 14a and 14b, respectively. 
The deformation causes dilatation and a corresponding out-of-plane stress, even though the 
deformation is nominally isochoric in the regions between the applied velocities, and all nodal 
velocities in the other two directions are set to zero. There are significant stresses in the vertical 
and horizontal directions, but these are primarily due to rotation of the stress tensor, and the 
magnitude is not attributed to kinematic coupling. The out-of-plane stress in figure 14a should be 
unaffected by the stress tensor rotation. The dilatation and the out-of-plane stress are not 
consistent with the boundary conditions and material constitutive response. 

It should be noted that the stress and dilatation are less in the uniform grid regions (but not 
adjacent to the edges where boundary effects play a role). This supports the assertion that these 
artifacts are created by the graded grid and would not occur if the grid were uniform. 
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Figure 13. Horizontal velocity field (m/s) creating two 
simple shear regions. The solid color regions 
at the top, center and bottom are nodes where 
the velocity is imposed. 
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Figure 14. The (a) out of plane stress (MPa) and (b) dilatation resulting from imposed simple shear boundary 
conditions on the graded mesh. 

3.2.5 Impact with a Spherical Projectile 

The sphere impact calculation of section 3.1.4 is repeated with ductile and brittle state-based 
materials. A 30-intn square target, 20 mm deep with a 0.5-tnm grid spacing, was impacted by a 
dense, rigid 10-mm diameter sphere traveling at 200 mm/ms. The target had a density of 
0.0022 g/mnr, a Young’s modulus of 70 GPa, a Poisson’s ratio of 0.28, and a yield strength of 
1300 MPa. The failure strain was set to 100 for the ductile material and to 0.01 for the brittle 
material. 
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The dilatation for the ductile target is shown in figure 15 at approximately 0.034 ms after impact. 
At this time, the penetration depth is over 6 mm—deeper than the penetrator radius. Due to the 
high failure strain, the damage is insignificant, on the order of l.e-6. There are two important 
observations from this image. First is that the surface heaving does not appear consistent with the 
volume of material displaced by the impactor (the vertical displacement on the sides is nearly 
zero). Given the low surface displacement, considerable densification would be expected, as in 
figure 7b. However, the levels indicated are much lower. The second observation is that the 
node locations and dilatation values under the penetrator are in disarray; some dilatation values 
even indicate expansion. Together, these observations suggest a breakdown in the ductile 
calculation at these very large deformations. 

The dilatation plots for the sphere impacting the brittle material at approximately 0.034 ms are 
shown in figure 16. The dilatation value is set to zero upon material failure. As with the ductile 
calculation, there appears to be insufficient surface displacement to accommodate the volume of 
the impactor. The nodes appear significantly compressed beneath the impactor and are likely 
constrained by the short-range interactions, since the specified constitutive relations would not 
apply for failed material. The region immediately below the sphere is punched vertically, and 
there are large velocity discontinuities associated with the fractures on either slide plug. 

1 -5.D00E-04 
-3.25GE-0 3 

-6.319E-0 3 
-9.381E-0 3 
-1.244E-0 2 
-1.551E-0 2 
I -1.857E-0 2 
-2. 163E-0 2 

1 -2.469E-0 2 
-2.776E-0 2 
-3.082E-0 2 
-3.388E-0 2 
-3■694E-0 2 
-4.D01E-02 
-4.307E-0 2 
-4.813E-02 
-4.950 E-Q2 



Figure 15. Dilatation plot for a dense, rigid sphere impacting a ductile target. The disordered 
nodes and color levels combined with an apparent densification suggest a 
breakdown of the calculation. The initial configuration is shown in figure 7a. 
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Figure 16. Dilatation plot for a dense, rigid sphere impacting a brittle target. The tight 
node packing suggests that the hydrodynamic response is being controlled 
by the short-range interaction force rather than the constitutive response. 


4. Kinematics of the State-based Formulation 


The results from EMU show volume-shear coupling when using both a shifted grid (figure 4) 
and a spatially graded grid (figure 12). Such couplings are inconsistent with the governing 
equations. It is instructive to investigate the behaviors to determine whether they result from the 
kinematics, the momentum balance, or some numerical aspect. The kinematic description is 
analyzed in this section through both analytic solutions and simplified numerical results. 

4.1 Analytic Solutions 

4.1.1 Normal impact with a Shifted Grid 

The volume-shear coupling in a shifted grid is examined by considering the 2D sheared mesh in 
figure 17. The vertical spacing of the horizontal lines of nodes is one unit, and the nodes are one 
unit apart along the horizontal lines. Each successive row of nodes is shifted by 0.2 units to the 
right from the row below it. A horizon of 2.3 units is used to keep the problem easily tractable. 
To simulate the initial stages of longitudinal wave propagation, a velocity is assumed on the 
bottom row of nodes and the remainder of the nodes is still at rest. 
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Figure 17. Shifted 2D grid used to investigate volume-shear coupling. 


The i, values for the labeled nodes are calculated as 

= (-2,0); ^ = (-1,0); ^ = (1.0); £, = (2.0); ? 4 = (-1.8,1); 

? s = (-0.8,1); ? 6 = (0.2,1); = (1.2,1); & = (-0.6,2); & = (0.4,2). 1 J 


These are combined to create the K tensor, as in equation 4, assuming that the node volumes are 
one and the weighting function is one within the horizon. The components of K and K 1 are 

9 



y (* „ ® f„) 



O' 


15.88 -1.61 

and 

Kf 1 - 1 

12 

1.6 ] 

L—1.6 12 J 


11 188 

1-1.6 

15.88^ 


(14) 


n=0 


Assume a vertical velocity on the first row of V; let v n be the velocity vector of the nodes and let 
v x be the velocity vector of the focus node. The Cartesian components of the velocity gradient 
are 
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The non-zero 21 component of the velocity gradient indicates a shear strain rate, which creates a 
shear stress. The 22 component of the velocity gradient would be equal to the velocity divided by 
the grid spacing in a local algorithm, but it is low by almost 34%. The nonlocal nature of the 
formulation encompasses additional nodes and spreads the displacement over a wider area. It is 
interesting to note that if nodes 8 and 9 are omitted from the calculations, the shear is zero and 
the error in the 22 component of the velocity gradient is also zero. Removing these points 
removes much of the nonlocal influence for this configuration, and it resembles a local model. 


With the exception of a regular grid and a normal velocity, it is not obvious how to discern 
which node arrangements will produce a shear deformation and which will not. In the shifted 
grid example, the shear deformation is dependent on the magnitude of the shift and the horizon. 
To illustrate this point, consider the horizon of 2.3 units in the previous example. If the nodal 
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shift per row is less than 0.06789 of the grid spacing, the shear will be zero. A horizon of 2.5 
times the grid spacing allows for a grid shift per row of 0.25 of the grid spacing without shear 
error. The reason for this is the number and symmetry of nodes entering or leaving the horizon. 


4.1.2 Point-wise Solution on Irregular Grid 


The effect of a node shift perpendicular to the wave can potentially be offset by shifts on other 
nodes. A mesh that appears severely distorted may not necessarily produce a shear deformation 
at every point. For example, consider the grid in figure 18. 
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Figure 18. Grid with single point shifted to create a zero shear error at a specific node. 


The mesh in figure 18 has an origin at (0,0) and a horizon of 2.2 units; all 13 digits are 
significant for this example. A boundary condition forces nodes on and below the x-axis to have 
constant velocity in the y-direction. Flere the nodal volume, V n , is calculated as 


V n = On - *o) 2 + On “ To Y ■ 


(16) 


The shape tensor is calculated as 
12 


K = 'Y d tfn®tn)*V n = 


14.7589408826158 0.717704414172386 
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and the inverse is 
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Letting V be the v-direction velocity of the nodes below the x-axis yields a velocity gradient at 
the origin of 


L = 


V 


0 

7.43849426498855 * 10“ 16 


0.3598602 999 369 3 1j 


(19) 


This illustrates that the shear defonnation for this mesh at this point is zeroed by judicious 
placement of a single node. This adjustment, however, likely adds to the shear at adjacent nodes. 


Two points can be gleaned from this example. First, the only way to be certain of the kinematic 
coupling that an irregular grid will generate in peridynamics is to calculate it. There is no 
obvious pattern or rule of thumb. Second, artifacts introduced by an irregular grid can vary from 
point to point and need not be smooth. 


4.1.3 Horizon Effect and Discretization Error 


The discretization error in computing the defonnation gradient can be estimated for simple 
configurations. Suppose a 1-D string of nodes is stretched such that the new position of the nodes 
is governed by the quadratic relation 

x = X + X 2 a . (20) 


Here a is a constant, x is the new position vector, and X is the original position vector. The ID 
deformation gradient is equal to the derivative of the displaced coordinate with respect to the 
initial coordinate. For this problem, the deformation gradient at any point equals 

dx 

— = 1 + 2 aX . (21 

ax 

According to the above equation, the defonnation gradient at the origin (X = x = 0) should be 
equal to 1. Assuming a horizon of 3, the nonlocal shape tensor is 28 and the deformation vector 
state is 28. This gives a deformation tensor of exactly 1 at the origin. In fact, any horizon used 
with this equation and a uniform grid will give the exact result at all nodes. 


Next, consider a higher order displacement where the new positions are governed by the cubic 
equation 

x = X + X 3 a . (22) 


For this displacement field, the deformation gradient at any point should be equal to 

dx 


dX 


= 1 + 3aX 2 


(23) 


The deformation gradient at the origin (X — x — 0) should again be equal to 1. With a horizon 
of 3, the nonlocal shape tensor is 28, but the deformation vector state is now 28 + 196a; this 
gives a deformation gradient of 1 + 7 ad 2 , where d is the node spacing. For any non-zero value 
of a, the deformation gradient has an error. This result is compared to the same problem with a 
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different horizon. Assuming a horizon of 2.9 times the grid spacing, the deformation gradient is 
1 + 3Aad 2 , and for a horizon of 1.1, the deformation gradient is 1 + ad 2 . This indicates that the 
error in the defonnation gradient is dependent on the chosen horizon as well as the grid spacing. 
It is not difficult to show that the error for the cubic displacement field is 

Y H - n 4 

Error — ad 2 _,^~ 1 —- , (24) 

A n =i n 

where H is the number of nodes included in the horizon to either side of the node of interest. 
Following this trend, the error is reduced as the horizon is decreased. When the horizon is 
reduced to a single node spacing, a local model rather than a nonlocal model, the ID 
peridynamics solution resembles the finite difference method. 


4.2 Numerical Evaluation 

The analytic results illustrate the effect of nodal placement on the calculation of the velocity 
gradient using the state-based formulation. These solutions are for individual points and do not 
give an overall picture of the distribution of the volume-shear coupling. To provide a larger scale 
evaluation, a stand-alone, 2D code was written to calculate the deformation gradient using the 
relations given by Foster, et al. (7). The initial and displaced coordinates are specified for all 
points, and the only aspect of the state-based peridynamics solution being evaluated is the 
calculation of the deformation gradient. 


4.2.1 Volume-shear Coupling 


The stand-alone calculations were used to investigate the volume-shear coupling for several 
mesh configurations: sinusoidal, graded, and rotated (figure 19). In all cases, the same 
displacement field is applied. The initial and final x-coordinates are identical. The displaced y- 
coordinates of the upper portion of the grid are the same as the initial coordinates, so the 
deformation gradient is the identity tensor in the upper region. The coordinates in the lower 
portion of the grid are shifted up 10% of the average grid spacing, and the deformation gradient 
for this rigidly displaced region is also the identity. The deformation gradient is nonuniform in 
the transition region. The transition follows a hyperbolic tangent function of the form 


y — Y + (0.1 d) - 


1 — tank 


M 


Y - 10d\y 

d jj. 


(25) 


where Y is the initial coordinate, y is the displaced coordinate, and d is the grid spacing. The grid 
spacing was 0.3 mm in the v-direction, and the mesh height was approximately 16 mm. 
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Figure 19. Grid configurations illustrating volume-shear coupling, (a), (c) and (e) show the volume change 

imposed on the grid. Plots (b), (d) and (f) illustrate the corresponding shear induced by the volume 
change. 


The results are shown in figure 19. The left column shows the grid and the volume change (Det F 
- 1). The corresponding plots at the right show the shear induced by this volume change. The 
shear induced in the sinusoidal mesh, shown in figures 19a and 19b, is anticipated because EMU 
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had already shown that a sheared mesh will induce the coupling. The coupling for the graded 
mesh, shown in figures 19c and 19d, was also observed in EMU, but figure 19d shows additional 
localized high shear spots that correspond to points dropping into and out of the horizon. The last 
results in figure 19 are for a rotated mesh. Here the sharp gradient induces pattern of positive and 
negative shear. 

These grids were also run as a wave propagation calculation in EMU. Some degree of volume- 
shear coupling was observed in all. The shear shown by EMU for the sinusoidal mesh 
strengthens and weakens as the mesh curvature changes. Volume-shear coupling is observed at 
early times for the rotated mesh, but as the wave front spreads over several grid spacings, the 
shear fades into the noise. The graded mesh calculations in EMU were presented earlier. 

It is important to remark that a linear displacement field (affine deformation) gives no shear for 
any of the grids shown in figure 19. This is consistent with the proof accompanying equation 13 
in Foster, et al. (7). In addition, the shear contribution is zero if the nodes are set on a regular 
orthogonal grid aligned with the wave propagation. 

An additional feature to note in figure 19 is that the shear deformation is highest at the 
boundaries, where the nodes are asymmetric within the horizon. 

4.2.2 Horizon Effects in a Collapsing Cylinder 

The second assessment of the peridynamic state-based defonnation gradient calculation was on a 
regular orthogonal grid, since most of the grid effects discussed thus far are greatly diminished 
for a regular orthogonal grid. The closed fonn analysis in section 4.1.3 suggests sensitivity of the 
solution to the displacement gradient and to how points are included in the horizon, so these 
aspects are investigated numerically. 

The problem considered is the isochoric, radial collapse of a cylinder. The analytic solution is 
simple, and the inhomogeneous deformation samples the grid at every angle. The cylinder has a 
20-mm outside radius and a 10-mm inside radius (A), and it is defonned modestly so that the 
inside radius is displaced to a radius of 9 mm (a). The regular grid spacing chosen was 0.3 mm. 
The deformed coordinates of all nodes were prescribed consistent with an isochoric radial 
collapse. The coordinates remain on a radial line with original radius R and the defonned radius 
given by 

r = Vfi 2 ~A 2 +a 2 , (26) 

which is based on area calculations and is valid to arbitrarily large strains. While only a 1/4 
segment of the cylinder is shown, a larger portion of the cylinder was analyzed to eliminate edge 
effects on the x and y symmetry axes. 

Results are shown in figure 20 for progressively larger horizons that envelop an additional node 
in rows and columns adjacent to the evaluated node. The horizons are 1.015, 1.5, 2.015, 2.3, 
3.015, and 3.2 times the grid spacing. The quantity plotted is the volume change (Det F - 1), 
which should be zero for the given isochoric deformation. The existence of discretization error is 
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not unexpected, although the magnitude is of some concern. The pressure error within the 
cylinder would be a significant fraction of the flow stress for a typical metal. 



Figure 20. Calculated volume change resulting from nodal displacements corresponding to an isochoric 

deformation field. The horizon for the plots is: (a) 1.015, (b) 1.5, (c) 2.015, (d) 2.3, (e) 3.015, and 
(f) 3.2 times the grid spacing. The color levels are the same for each plot. 
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What is interesting in figure 20 is how the patterns and magnitude of the error shift as the 
horizon changes to include more points. The changes occur in steps, so a point moving into or 
out of the horizon would cause a step change in pressure. Also significant is the large volume 
change at the inner and outer radii. Here, the volume strain can exceed half a percent; it occurs 
abruptly in transition from the interior to the boundary regions; and multiple transitions occur in 
the boundary region. 


5. Discussion 


One of the frequently cited features of peridynamics is that it replaces partial differential 
equations of conventional continuum theory with integral equations ( 1, 8, 9). However, when the 
peridynamic and FE strain calculations are viewed side by side at comparable stages of their 
formulations, this distinction is not apparent. Both involve integrals over volumes. The cited 
distinction of integral versus differential equations is a result of comparing different aspects of 
the formulations. 

The integrands for the two methods do differ. The conventional continuum methods use analytic 
derivatives, and peridynamics uses a difference approximation to the derivative. At a small 
enough point spacing, the difference approximations to gradients can also resemble singularities. 
As an example, this construct can be taken to the smallest sensible level—the atomic spacing. If 
the material points are taken to correspond to atomic positions with initial spacing on the order of 
2.5 x 10 4 microns, a 100-micron-wide crack would have a stretch of 4 x 10 5 ; this is not the 
infinite stretch associated with a crack in a traditional continuum framework, but it is still large 
enough to dominate the behavior of the integral. 

This distinction in gradient representation becomes academic when the equations are discretized 
through nodal relations. The deformation gradient and velocity gradient are formed from sums 
and differences among nodal coordinates and velocities. The manner in which the nodal values 
are combined can vary considerably for the two methods, and, as a nonlocal method, 
peridynamics includes significantly more nodes in the strain calculations. Consequently, the 
results can be quite different. 

The nodal force calculations in peridynamics are significantly different from FE nodal force 
calculations. The momentum balance in peridynamics uses a difference between stress vector 
states at neighboring nodes. This resembles the strong fonn of the momentum equations, 
whereas FE typically uses a weak fonn. 

Another cited attribute of peridynamics is that the fracture occurs naturally on bonds connecting 
nodes rather than through additional equations, stress intensity factors, and prescribed crack 
paths that must be imposed in conventional continuum fonnulations (1 ). This comparison is to 
special techniques used in simulations where detailed crack stress fields are resolved, rather than 
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the more broadly used FE failure methods. The failure criteria in peridynamics are quite similar 
to simple FE failure models used in large-scale hydrocodes. The peridynamics micro models use 
bond stretch to trigger failure, and the state-based models use a critical strain to initiate failure. A 
strain to failure is one of the simpler failure criteria applied in FE analyses. 

Material failure in peridynamics is progressive and should exhibit less grid sensitivity than the 
abrupt failure that is often used in finite elements. Detailed comparison of peridynamics fracture 
and FE failure using nonlocal constitutive relations or progressive stress release algorithms is 
needed. 

The post-failure tensile response of peridynamics allows unrestrained separation of nodes on 
either side of the fracture surface. Lagrangian FE allows a similar unhindered separation, 
although element deletion may be necessary to prevent numerical issues when there are large 
displacements between the two surfaces. 

Material densification in peridynamics following compressive failure is controlled by the short- 
range interaction forces. These forces result from separate equations, not the constitutive 
relations. It is possible to have unhindered compression between the time of fracture and the 
activation of the short-range force. This unhindered compression is primarily responsible for the 
unrealistic densification seen in the impact simulations. In conventional Lagrangian, FE the 
pressure-volume relation within an element is still enforced following failure if the pressure is 
positive. In these cases, conversion of failed material to smooth particle hydrodynamics (SPH) 
particles (12) can be used to provide continuous repulsive forces at large defonnations. It should 
be possible to modify peridynamics to use the equation of state so that a separate short-range 
interaction force is not necessary. It might also be possible to modify the short-range forces in 
peridynamics to achieve a continuous compression response. 

As links between nodes are removed due to failure, the distribution of nodes within the horizon 
contributing to the solution will change. It was demonstrated previously, e.g., figure 20, that such 
changes to the node distribution can induce step changes in nodal strain fields, even if the node is 
not immediately adjacent to the fracture. The volume and/or shear defonnation can experience 
jump discontinuities when nodes are added to or removed from the horizon. The discontinuous 
nature of the mesh effects makes it difficult to identify and quantify errors. The magnitude and 
effect of such jumps on the details of the stress and strain fields during failure should be 
examined to determine whether or not the crack propagation solution is affected significantly. 

An additional aspect of failure to consider is the fracture surface. The nonlocal fonnulation 
introduces boundary effects that will be present on these newly created free surfaces. Figures 8, 
14, 19, and 20 all show effects associated with surfaces. In figure 8, the surface nodes and the 
nodes 3 rows from the surface have nearly identical displacement histories for a wave transit 
problem. This is not expected from the governing equations. In figures 14 and 19, the edges have 
larger volume-shear coupling than the interior, and in figure 20, an artificial volume change of 
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more than 0.5% was calculated at the boundary for the isochoric displacement field. There is no 
reason to expect that these surface effects would not also be present at crack surfaces. 

In light of the combination of nonlocal effects at free surfaces and effects from li nk s broken 
within the horizon, the accuracy of the stress and strain solutions near a crack should be 
examined in detail. The crack tip displacement, strain, and stress fields should be compared to 
analytic solutions to determine the extent to which the nonlocal formulation and numerical 
discretization impact the results. This is beyond the scope of the current study. 

Difficulties associated with boundary conditions in nonlocal methods may also be important for 
other surfaces and interfaces. If the strain or stress state is perturbed significantly at a surface, 
evaluation of the failure criteria at the surface could also be affected. Alternatively, if two 
surfaces impact and the contact is only made at the surfaces nodes, the impact conditions may 
not be imposed properly. In order to create the expected surface behaviors, it may be necessary 
to refine the grid sufficiently to reduce the spatial extent of the perturbations or to develop 
special boundary conditions consistent with the numerical treatment, as is done in finite 
difference implementations (13), certain hydro codes (e.g., CTH [14)), and other nonlocal 
methods, such as Element Free Galerkin (15). 

Many results were presented illustrating grid-dependent volume-shear coupling for irregular and 
non-orthogonal grids. The effects result from a combination of the nonlocal method and 
asymmetries of the node distribution within the horizon. Such effects are not consistent with the 
governing equations and will introduce unanticipated shear and pressure into the results. One 
possible way around these effects may be to only use regularly spaced orthogonal grids. These 
appear to avoid the most severe volume-shear coupling. However, there is still some coupling for 
inhomogeneous deformation fields, as illustrated by the collapsing cylinder; and it is not possible 
to create conformal boundaries on arbitrary geometries using regular grids. Further detailed 
evaluations are necessary to determine conditions under which peridynamics will produce 
accurate general defonnation results, including strain and stress states at boundaries and 
interfaces. 

Dissipation is a significant aspect of the peridynamic solution in the context of wave 
propagation. Stress waves often set the loading for high-rate fracture simulations, so wave 
propagation can be quite important. The nonlocal nature of the strain and momentum 
calculations averages the fields within the horizon. The horizons from adjacent nodes have 
significant overlap and this overlap can propagate infonnation to neighboring nodes. This was 
illustrated by the effect of horizon on rise time seen in comparing figures 2a and 3a. Dispersion 
is also evidenced by the variation in wave propagation speed, where the lower particle velocity 
moves faster through the material than the higher particle velocity. This velocity difference leads 
to further spreading of the wave front. Additional dissipation comes from application of drag 
forces to quiet noise in the solution. The drag algorithm ties all nodes within a family, thereby 


29 



suppressing peak values and further diminishing gradients, as shown in figure 9. It has a much 
broader effect than control of zero-energy modes. 

The long-range interactions in peridynamics necessitate many equations for each node. A typical 
concern for techniques where the equations greatly outnumber the degrees of freedom is that 
solutions will be over-constrained. This will occur if there is inconsistency among the constraints 
arising either from the underlying equations or the numerical approximation. The possibility that 
over-constraint is responsible for the unexpected velocity patterns in figure 8a and the node 
locations in figure 15 should be investigated. 

A final discussion topic concerns solution accuracy in relation to the number of points in the 
horizon. The relative position vector, is numerically larger for nodes closer to the horizon than 
for nodes near the center of the integration region. Multiplying the differences in field quantities 
by the relative position vector effectively gives a greater weight to values from more distant 
nodes. Concurrently, including the distant nodes increases the magnitude of the shape tensor. 
Dividing by the larger shape tensor also reduces the influence of the nodes adjacent to the node 
of interest. The combined result is that accuracy can be diminished significantly as the horizon is 
increased, as illustrated in section 4.1.3. Thus, it would seem that the most accurate solutions 
would be obtained by only including adjacent nodes. This would result in a local model, and the 
discretized peridynamics equations would resemble the well-known finite difference relations. 
Alternatively, the weighting function could be employed to ensure that more distant nodes have 
progressively less impact, as in element free Galerkin methods (15). This is another area for 
future study. 


6. Conclusions and Recommendations 


A side-by-side comparison of defonnation gradient calculations for peridynamics and finite 
elements suggests that the oft-cited distinction of integral equations versus differential equations 
is likely an impression from comparing different aspects of the formulations—both methods use 
integrals over volumes, but the integrands use different approximations. The greatest distinction 
is that peridynamics is a nonlocal method, whereas traditional finite elements is a local method. 

A survey of the peridynamics solutions for pre- and post-fracture mechanical defonnation was 
undertaken using a combination of analytic methods and the EMU and KRAKEN codes from 
Sandia. The study revealed anomalous behaviors related to the numerical grid, the nonlocal 
formulation, and the particular model implementations. 

The nonlocal method and grid effects can produce artifacts that are not consistent with the 
governing equations: 
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• Volumetric and shear defonnations are coupled through the kinematics if the node 
distribution within the horizon is irregular. A volume change can induce shear deformation, 
and a shear defonnation can induce a volume change. 

• The solution can deviate considerably from expectation on boundaries because nonlocal 
boundaries have a finite thickness of several node spacings. Fracture surfaces are also 
boundaries and may be impacted by nonlocal effects. 

• There is high dispersion of field quantities, most notably, significant spreading of wave 
solutions. 

There are also features related to the particular implementation and material models used. These 
result from implementation choices and should be relatively straight-forward to improve. 

• When the micro-elastic and micro-plastic models are stressed beyond the yield point the 
results are not consistent with the physical behavior of materials. The pressure-volume 
response also follows the plasticity relation and the effective bulk modulus is softened. 

• Compression response of state-based materials at large defonnation and following fracture 
is governed by the short-range interaction force rather than the equation of state. 

• The drag algorithm controlling zero-energy modes is applied broadly as a smoothing 
function rather than explicitly constraining only hourglass modes. Specific zero-energy 
modes are not identified. 

• Solution errors can grow significantly as the integration horizon is increased. 

Given the current modeling capabilities of EMU and KRAKEN, potential application for the 
peridynamics codes are brittle tensile fracture for the micro-elastic models and more general 
tensile fracture for the state-based models. Even for these applications, analysts should 
understand how the problem set-up and results may differ from more conventional approaches, 
and ensure that the behaviors are consistent with expectations, particularly at boundaries. 
Algorithm modification and extensive, detailed validation of defonnation and stress fields would 
be needed for applications involving fracture and compression, particularly for large strains 
associated with high-speed impact simulations. 

As with any new computational tool, the analyst should run numerous detailed comparisons with 
existing solutions, analytic methods, and data on the class of problems of interest to verify that 
the results are conect and that the new method has quantifiable advantages. The benefits of 
potentially less mesh sensitive fracture have to be balanced against grid sensitivities, boundary 
effects, and limited constitutive behaviors. 
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